#include <iostream>
#include <fstream>
#include <vector>
#include <memory>
#include <string>
#include <cassert>
#include <direct.h>

using namespace std;

string GetInputDirectory(string fileName);
string GetOutputDirectory(string fileName);
double CalculateRe(double rhoG, double DiaP, double slipVel, double muG);
double CalculateCd(double Re);
double CalculateTao(double rhoP, double DiaP, double muG, double Cd, double Re);

int main()
{
	//////////////////////////////////////////////////////
	// Section 1: load in data
	//////////////////////////////////////////////////////
	string inputFileName = "inputData.txt";
	string inputFileDirect = GetInputDirectory(inputFileName);
	ifstream inputFile(inputFileDirect);
	double pX0 = 0.0;
	double dt, pV0, pRho, pDia, velG, muG, rhoG;
	assert(inputFile.is_open());

	string tempString;
	inputFile >> tempString >> dt;
	inputFile >> tempString >> pV0;
	//inputFile >> tempString >> Cd;
	inputFile >> tempString >> pRho;
	inputFile >> tempString >> pDia;
	inputFile >> tempString >> velG;
	inputFile >> tempString >> muG;
	inputFile >> tempString >> rhoG;
	inputFile.close();

	/////////////////////////////////////////////////
	 //Section 2: 
	//Calculate the velocity and distance at each time
   //////////////////////////////////////////////////
	int nTsteps = 10;
	auto tAll = make_shared<vector<double>>(nTsteps + 1, 0.0);
	auto pVAll = make_shared<vector<double>>(nTsteps + 1, 0.0);
	auto pXAll = make_shared<vector<double>>(nTsteps + 1, 0.0);
	(*pVAll)[0] = pV0; (*pXAll)[0] = pX0;
	for (int i = 0; i < nTsteps; i++)
	{
		int i_store = i + 1;
		(*tAll)[i_store] = (*tAll)[i] + dt;
		double pReTemp = CalculateRe(rhoG, pDia, abs(velG - (*pVAll)[i]), muG);
		double CdTemp = CalculateCd(pReTemp);
		double pTaoTemp = CalculateTao(pRho, pDia, muG, CdTemp, pReTemp);
		//Update particle velocity
		(*pVAll)[i_store] = (*pVAll)[i] + dt * (velG - (*pVAll)[i]) / pTaoTemp;
		(*pXAll)[i_store] = (*pXAll)[i] + dt * (*pVAll)[i];
	}

	/////////////////////////////////////////////////
	 //Section 3: 
	//Write out data to file
   //////////////////////////////////////////////////
	string outputFileName = "outputData.txt";
	string outputFileDirect = GetOutputDirectory(outputFileName);
	ofstream outputFile(outputFileDirect);
	assert(outputFile.is_open());
	for (int i = 0; i < nTsteps + 1; i++)
	{
		if (i == 0)
		{
			outputFile << "time (s) " << "Vp (m/s) " << "X (m) " << endl;
		}
		outputFile << (*tAll)[i] << " " << (*pVAll)[i] << " " << (*pXAll)[i] << endl;
	}
	outputFile.close();

	//////////////////////////////////////////////////////////
	//clear pointer
	tAll.reset(); pVAll.reset(); pXAll.reset();

	return 1;
}

string GetInputDirectory(string fileName)
{
	char homeDirect[1024];
	_getcwd(homeDirect, 1024);
	string localDirect = (string)(homeDirect)+(string)("\\LocalSettings\\") + fileName;
	return localDirect;
}
string GetOutputDirectory(string fileName)
{
	char homeDirect[1024];
	_getcwd(homeDirect, 1024);
	string localDirect = (string)(homeDirect)+(string)("\\LocalResults\\") + fileName;
	return localDirect;
}
double CalculateRe(double rhoG, double DiaP, double slipVel, double muG)
{
	return rhoG * DiaP * slipVel / muG;
}
double CalculateCd(double Re)
{
	// empirical spherical model by Mori and Alexander
	double Cd = 0.0;
	if (Re < 0.1)
	{
		Cd = 24.0 / Re;
	}
	else if (Re >= 0.1 && Re < 1.0)
	{
		Cd = 22.73 / Re + 0.0903 / (Re * Re) + 3.69;
	}
	else if (Re >= 1.0 && Re < 10.0)
	{
		Cd = 29.1667 / Re - 3.8889 / (Re * Re) + 1.222;
	}
	else if (Re >= 10.0 && Re < 100.0)
	{
		Cd = 46.5 / Re - 116.67 / (Re * Re) + 0.6167;
	}
	else if (Re >= 100.0 && Re < 1000.0)
	{
		Cd = 98.33 / Re - 2778.0 / (Re * Re) + 0.3644;
	}
	else if (Re >= 1000.0 && Re < 5000.0)
	{
		Cd = 148.62 / Re - (4.75e4) / (Re * Re) + 0.357;
	}
	else if (Re >= 5000.0 && Re < 10000.0)
	{
		Cd = -490.546 / Re + (57.87e4) / (Re * Re) + 0.46;
	}
	else if (Re >= 10000.0 && Re < 50000.0)
	{
		Cd = -1662.5 / Re + (5.4167e6) / (Re * Re) + 0.5191;
	}
	else
	{
		assert(false && "Cd value is out of Range!");
	}
	return Cd;
}

double CalculateTao(double rhoP, double DiaP, double muG, double Cd, double Re)
{
	return rhoP * pow(DiaP, 2) * 24.0 / (18.0 * muG * Cd * Re);
}